Non-linear evolution of f{R) cosmologies II: power spectrum 
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We carry out a suite of cosmological simulations of modified action f{R) models where cosmic 
acceleration arises from an alteration of gravity instead of dark energy. These models introduce 
an extra scalar degree of freedom which enhances the force of gravity below the inverse mass or 
Compton scale of the scalar. The simulations exhibit the so-called chameleon mechanism, necessary 
for satisfying local constraints on gravity, where this scale depends on environment, in particular the 
depth of the local gravitational potential. We find that the chameleon mechanism can substantially 
suppress the enhancement of power spectrum in the non-linear regime if the background field value 
is comparable to or smaller than the depth of the gravitational potentials of typical structures. 
Nonetheless power spectrum enhancements at intermediate scales remain at a measurable level for 
models even when the expansion history is indistinguishable from a cosmological constant, cold 
dark matter model. Simple scaling relations that take the linear power spectrum into a non-linear 
spectrum fail to capture the modifications of f{R) due to the change in collapsed structures, the 
chameleon mechanism, and the time evolution of the modifications. 
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I. INTRODUCTION 

Cosmic acceleration can arise from either an exotic 
form of energy with negative pressure or a modification 
to gravity in the infrared. Self-consistent models for the 
latter are highly constrained by the stringent tests of 
gravity in the solar system. Additional propagating de- 
grees of freedom must be suppressed by non-linearities 
in their equations of motion [l|, Q, 0] in a stable manner 
05 S 01 • These suppression mechanisms manifest them- 
selves with the formation of non-linear structure in the 
Universe 0, H, @|- Understanding the physical content, 
phenomenology and even the basic viability of such mod- 
els thus requires cosmological simulations. 

One possibility that has received much recent atten- 
tion is the so-called f{R) class of models (see [10| and 
references therein). These models generate acceleration 
through a replacement of the Einstein-Hilbert actio n by 
a function of the Ricci or curvature scalar R [ill , [l2l , [l3( . 
They also introduce an extra propagating scalar degree 
of freedom that acts as an effective fifth force on all forms 
of matter [ij, [lH . The range of the force depends non- 
linearly on the local curvature and can be made to be- 
come infinitesimal at high curvature. With an appropri- 
ate choice of the function f{R)^ deep potential regions 
can trap the field at high curvature leading to a non- 
linear "chameleon mechanism" |16| t hat suppresses local 
deviations from ordinary gravity [13, [l^, ^^ • 

Whether or not solar system tests of gravity are satis- 
fied in an f{R) model then depends on the depth of the 
gravitational potential including the astrophysical and 
cosmological structure surrounding it [19]. Likewise ob- 
servable deviations from ordinary gravity for upcoming 
dark energy probes such as weak lensing, galaxy cluster- 
ing and clusters of galaxies depend on the whole history 
of non-linear structure formation. 

Previous cosmological simulations {e.g. [2Q|, \2u. [221, [231, 



[2^ ) have focused on modifications of the force law with a 
fixed and density independent range. Such modifications 
alone are incapable of satisfying local tests of gravity. 

In a companion paper [25], the numerical methodology 
for solving the non-linear field equation of f{R) gravity 
was established. In this second paper of the series, we 
apply this methodology and carry out a suite of cosmo- 
logical simulations of f{R) models that are chosen to 
expose the impact of the chameleon mechanism on the 
power spectrum of the matter and the lensing potential. 

We begin in ^TT] with a review of non-linear gravita- 
tional dynamics in f{R) models and proceed to the sim- 
ulation results in ^III[ We discuss these results in ^IV[ 



II. f{R) DYNAMICS 

A. Basic equations 

The class of modifications we consider generalizes the 
Einstein-Hilbert action to include an arbitrary function 
f{R) of the scalar curvature R 



S 



d x\ 



R + f{R) 
167rG 



(1) 



Here Lm is the Lagrangian of the ordinary matter which 
remains minimally coupled. Setting f{R) = recovers 
general relativity (GR) without a cosmological constant 
whereas setting f{R) = —WnGpA = const, recovers it 
with a cosmological constant. Here and throughout c = 
h=l. 

Variation of Eq. ^ with respect to the metric yields 
the modified Einstein equations 



Gaf3 + Fa(3 = SirGTap , 



(2) 



-Fa/3 = fuRaff ~ ( 2 ~ ^^^ I ^"'^ ~ ^oypfR , 



where the field, 



fR 



dfjR) 
dR ' 



(3) 



plays the role of a propagating extra scalar degree of 
freedom. In particular the trace of the modified Einstein 
equation ([2]) yields the equation of motion for the field 



□/fi = 



dfR 



(4) 



with the effective potential defined by 



^^\[R-fr,R+2f- SirGip - 3p)] . (5) 



The effective potential fias an extremum at 
R-RfR + 2f = SnGip - 3p) , 
and its curvature is given by 






i/dR 



R 



(6) 



(7) 



This can be interpreted as the effective mass of the field 
Jr and defines the range of the force. For stability, the 
extremum should be a minimum and hence /i^ > [5]. 

Phenomenologically viable models typically must have 
very fiat f{R) functions such that |/i^| <C 1 at cosmolog- 
ical curvatures and larger. The model we simulate is in 
the class of f{R) functions proposed in [1^, 



/(i?)cx 



R 



AR 



(8) 



where A is a constant with dimensions of length squared. 
In the limit i? ^ 0, f{R) ^ as with GR with no cos- 
mological constant. For sufficiently high curvature that 
AR ^ 1, f{R) can be approximated as a constant, which 
drives the acceleration, plus a term that is inversely pro- 
portional to curvature. In this limit, we can approximate 
Eq. (O by 



f{R)^-WnGpA-fRO^. 



(9) 



where we have set the constant A to match some effective 
cosmological constant pA- Here we define ^o = R{z = 0) 
as the background curvature today and /rq = fniRo)- 

Taking \fRo\ <C 1, the background expansion follows 
the ACDM history with the same pA to leading order in 
|/i?o| uM' I^ particular the background curvature may 
be approximated as 



R^3H^[nm{l^zf^4nA] , 



(10) 



where Cli = S7TGpi{z = 0)/3Hq. We can also simplify 
the mass term in Eq. ^ p ^ {3dfR/dR)~^^'^ defining 



the comoving Compton wavelength or range of the field 
Ac as 



Ac 



M 



m 



R0\ 



.Rq 

R^' 



(11) 



Notice that the range of the interaction has a steep in- 
verse dependence on the local curvature R. 

We take the WMAP3 Q fiat background cosmology 
throughout: 1]^ = 0.24, Qa = 0.76, ^b = 0.04181, Hq = 
73 km/s/Mpc, and initial power in curvature fiuctuations 
As = (4.52 X 10"^)^ at A: = O.OSMpc-^ with a tilt of 
Us = 0.958. For ACDM these parameters give as = 0.76. 

For these values, the Compton wavelength for the 
background today is 



Aco ^ 3.2 



1/ 



RO 



10- 



Mpc. 



(12) 



Note that even for field amplitudes as low as |/i^o| ^ 10~^ 
where deviations from ACDM in the expansion history 
are comparably negligible, the modifications to gravity 
at the cosmological background density are order unity 
at astrophysically interesting scales of a Mpc [l9| . 

Furthermore, the field equation (|4]) can be simplified by 
neglecting the small /rR term and assuming that matter 
dominates over radiation, thus resulting in 



□/i? = J[i^-8^G(p^ + 4pA)], 



(13) 



Subtracting off the background values for the field, cur- 
vature and density yields 



BSfR = - [SR{fR) - SirGSprn] , 



(14) 



where SR = R- R, S/r = /r - fR{R), Sp^ = pm - Pm- 
Note that the field fiuctuation is defined by subtracting 
off the field evaluated at the background curvature and 
not the spatially averaged field value. The procedure 
eliminates the potential ambiguity of defining fiuctua- 
tions with a highly non-linear field equation. 

The minimum of the effective potential is at a curva- 
ture corresponding to the GR expectation SR = STrGSpm- 
If the field value achieves this minimum then the high 
density regions have high curvature and short Compton 
wavelengths which suppress the deviations from ordinary 
gravity. However, we shall see in TO Bl that whether the 
field value achieves this minimum at any given location 
depends on the depth of the gravitational potential. 

Finally, we work on scales much less than the horizon 
such that the quasi-static limit applies where time deriva- 
tives may be neglected compared with spatial derivatives. 
This corresponds to assuming that the field instanta- 
neously relaxes to its equilibrium value. More specifi- 
cally, relaxation must occur on a time scale that is short 
compared with the non-relativistic motion of particles. 
This approximation should be excellent for wavelengths 
that are not orders of magnitude larger than the Comp- 
ton scale below which field perturbations propagate near 



the speed of light. The field equation in comoving coor- 
dinates then becomes a non-linear Poisson-type equation 



V^^/fl = - [5R{fR) - SttGSp, 



(15) 



An explicit consistency test of this approximation in the 
cosmological context is given in [25[ . 

Since f{R) is a metric theory of gravity, particles move 
in the metric or gravitational potential in the same way 
as in general relativity. However the field acts as a source 
that distinguishes the two potentials in the metric 



ds' 



-{1^2^)dt^^a^{1^2^)dx^ 



(16) 



In the quasi-static limit, the modified Einstein equations 
([2]) imply that the Newtonian potential ^, whose gradi- 
ent is responsible for the motion of particles, is give n by 
a Poisson equation that is linear in 6pm and 5R [23, 



V^^ 



WnG 



a^Spm - -—SRifn) . 



(17) 



Equations (p!5|) and (p!7|) define a closed system for the 
gravitational potential given the density field. It is inter- 
esting to note that braneworld modified gravity models 
|28| obey a similar system of equations except that the 
effective potential involves field gradients ^]. 



B. Non-linear chameleon 

Before proceeding to the numerical solution of these 
equations, it is worthwhile to examine the qualitative 
aspects of the f{R) system of equations to expose po- 
tential observational consequences. In particular, these 
equations exhibit the chameleon mechanism under which 
modifications to gravity become environment dependent. 

The force law modifications are manifested by the ap- 
pearance of the second potential $ in Eq. (p!6|) . Grav- 
itational lensing and redshifts of photons depend on a 
combination of the two potentials 



,($-*) 



= -AnGa^Spn 



(18) 



Note that this relationship between the lensing poten- 
tial and the matter density is unaltered from the GR 
expectation. It is the Poisson equation for ^ that has 
an altered relationship to the matter density and gov- 
erns the motion of non-relativistic particles. Ordinary 
gravity is recovered if SR = STrGSpm^ and the equation 
for the gravitational potential reduces to the unmodified 
equation. 



V^^ = AirGa'^Spn 



and 



($-^) 



^ 



-^. 



(19) 



(20) 



Deviation from this relation locally are constrained by 
solar system tests of gravity at the level of |(<l> + ^)/^| < 
10~^ [2^. These deviations are related to changes in the 
field through the field equation ([15]), the Poisson equation 
([T7|) and the lensing potential ([T8|) 



V2($ + M/) = -v25/^ 



(21) 



If the force modifications are small, then |($ + ^)/^| <C 1 
for local contributions to the potential. 

This relationship gives a rule of thumb for the appear- 
ance of the chameleon effect. Consider an isolated spher- 
ically symmetric structure embedded in the background 
density. At the center of the object, the change in the 
field is related to the total depth of the potentials as 



A/,, = -(^ + ^). 



(22) 



as long as both remain finite at the center. 

To drive the range of the force modification to a scale 
much smaller than the background value, the field am- 
plitude Ifnl must be substantially below its background 
value l/^l at the center. For example, at the current 
epoch this requires 



|A/k|~|/koI~I* + *I«I*I 



(23) 



This is a necessary condition for the appearance of a 
chameleon. A chameleon also does not appear if most of 
the mass contributing to ^ is on scales above the Comp- 
ton scale of the background. In this case |<l> + ^| remains 
much smaller than ^ in this regime and the change in 
potential does not contribute to a substantial change in 
the field A/r. Thus the smooth potential contributed 
by very large scale structure does not enter into the 
chameleon suppression. 

If the local potential is not sufficiently deep, the field 
prefers to remain smooth at the background value and 
not track SR{fR) = SirGSpm- Energetically, the cost 
of high field gradients prevents the field from Mng at 
the local minimum of the effective potential [l9|. In 
this case the field remains near its background value 
and SR{fR) <C SirGSpm- Eq. ([ITJ) then implies that the 
strength of gravity is a factor of 4/3 greater than ordinary 
gravity. 

This consequence can alternatively be seen directly 
through Eq. (fT5|) . If the field fiuctuation is small then 
SR ~ {dR/dfR)Sfii and the field equation can be solved 
in Fourier space as a linear Poisson equation 






(24) 



where p = p{R) and defines the Compton scale in the 
background through Eq. ([TT]) . This solution combined 
with the Poisson equation (fT7|) gives 



e^{k) 



-^^^'V^F 



M^a^ 



/x^o^ 



a^Spmi^) , (25) 
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FIG. 1: The mean power spectra A^ = k^P{k)/2'n'^ of the cosmo- 
logical simulations without Jr modifications (ACDM or |/ko| = 0)- 
Upper panel: Average power spectrum of all simulations in com- 
parison to linear theory and the Smith et al. non-linear fit (solid 
line). Lower panel: Relative power spectrum offset from Smith et 
al. fit. The individual simulation box sizes averages are plotted in 
solid, long-dashed, and dashed lines. Results converge at approxi- 
mately the half-Nyquist wavenumbers (vertical lines). Points with 
errors show the mean and the 1-cr error bars computed using the 
bootstrap method, weighted by simulation box volume, out to the 
half-Nyquist wavenumbers (see text). 



which has an effective G that depends on wavelength such 
that G -^ 4/3G below the Compton wavelength in the 
background. Note that linearity in the density field is 
not assumed. Hence we will use this field linearization 
approximation to test the effects of the chameleon in the 
simulations. 

In summary order unity modifications to gravity are 
expected on scales smaller than the Compton scale of the 
background but away from the centers of deep gravita- 
tional potential wells of cosmological structure. In these 
regions, the chameleon mechanism suppresses the local 
Compton scale and hence the force modifications. 



III. f{R) SIMULATIONS 
A. Simulation description 

To solve the system of equations defined by the mod- 
ified Poisson equation (pT|) and the quasi-static Jr field 
equation (p!5|) in the context of cosmological structure 
formation, we employ the methodology described in [25^. 
Briefiy, the field equation for fji is solved on a regular 
grid using relaxation techniques and multigrid iteration 
[30, llil • The potential ^ is computed from the density 
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FIG. 2: Relative power spectrum enhancement over ACDM at 
a = 1 for the full fji simulation compared with the no-chameleon 
fpi simulations and linear theory. At high k, linear theory un- 
derestimates the absolute power for both ACDM and fji while 
overestimating the relative enhancement. Without the chameleon, 
power is sharply enhanced on scales smaller than the Compton 
scale in the background which increases with \fRo\ (see Eq. ([12])). 
For l/i^ol = 10~^ the chameleon strongly suppresses these enhance- 
ments at high k. For |/ho| = lO"'^, this suppression is nearly ab- 
sent except for a residual effect from the chameleon at high redshift. 
The 3 highest points in k have increased sampling errors since they 
utilize only the 64:h~^ Mpc boxes. Note that the linear prediction 
is the fractional enhancement in the linear power spectrum, that 

^^' (^/k, linear " ^ACDM, linear )/^ACDM, linear • 



and /r fields using the fast Fourier transform method. 
The dark matter particles are then moved according to 
the gradient of the computed potential, — V^, using a 
second order accurate leap-frog integrator. 

Since the Compton wavelength or range of the /r field 
in Eq. (pTj) shrinks with increasing curvature, modifica- 
tions to the force law and structure formation above any 
given comoving scale vanish at sufficiently high redshift. 
We can therefore treat the initial conditions for the sim- 
ulations in the same way as in a cosmology with ordi- 
nary gravity. All simulations have the cosmological pa- 
rameters given in the previous section which make them 
all compatible with the high redshift CMB observations 
from WMAP [32!] . To explore the modifications induced 
by the /r field, we simulate models with field strengths 
in the background of \fRo\ = 0, 10-^ 10-^ IQ-"^. Note 
that the \fRo\ = is exactly equivalent to ACDM. 

Specifically, the initial conditions for the simulations 
are created using Enzo [33], a publicly available cosmo- 
logical N-body + hydrodynamics code. Enzo uses the 
Zel'dovich approximation to displace particles on a uni- 
form grid according to a given initial power spectrum. 
We use the initial power spectra given by the transfer 
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FIG. 3: 2D slices through simulations at a = 1 for various fields: matter overdensity 6 = 6pm/pm (iGft); minimum gravitational potential 
^ along the line of sight (middle) and minimum field fn/fRo along the line of sight (right). Each slice has the two dimensions of full 
64/i~^ Mpc box and is projected across 16h~^ Mpc along the line of sight. Shading scale ranges from white to black between the extreme 
field values as follows: (-5.0 -^ 9.8) for ln(l + (5); (0.3 -^ -4.3) x 10-^ for ^; (1.05 -^ 0) for fn/fm. In the Jrq = \l^-^\ (top) simulation, 
the chameleon mechanism suppresses the field and hence the force deviations in the deep potential wells surrounding large overdensities. In 
the fjiQ = |10~^| (bottom) simulation, the fji remains stiff and a chameleon does not appear at the present. While the density and potential 
are slightly enhanced in this run due to the force modification, the appearance of the chameleon depends mainly on the background field 
value /ho- 



function of Eisenstein & Hu [34| and normalization fixed 
at high z. Note that the initial spectrum does not include 
the effects of baryon acoustic oscillations. The simula- 
tions are started at z = 49, and are integrated in time in 
steps of Aa = 0.002. 

In order to extend the dynamic range of the results, 
we run three simulations with box sizes Lbox = 256 h~^ 
Mpc, 128 h-^ Mpc, and 64 h'^ Mpc. Ah simulations 
are run with 512 grid cells in each direction and with 
A^p = 256^ particles. Thus, the formal spatial resolu- 
tions of the simulations are 0.5 h~^ Mpc, 0.25 h~^ Mpc, 
and 0.125 h~^ Mpc for the largest, middle, and small- 
est boxes, respectively. The corresponding mass resolu- 
tions are 2.76 x lO^^h'^ M©, 3.45 x lO^^h'^ Mq, and 
4.31 X lO^h-^ Mq. 

In the subsequent analysis, the particle Nyquist wave 
number of each simulation will play important roles in 
determining the range of trustworthy scales. As shown 
in [2J, [35[ , power spectrum of cosmological simulations 
start to show systematic > 10% deviations from Smith et 
al. 1361 and Peacock & Dodds |37| fits at wave numbers 



above half the particle Nyquist mode, defined as k^ = 
irNp /(2Lbox)- For our three simulation sizes, the half- 
Nyquist wave numbers are, in order of decreasing box 
size, 0.79/1 Mpc-\ l.bTh Mpc-\ and 3.Uh Mpc-^ 

Finally, to assess the impact of the chameleon mecha- 
nism on the power spectrum, we also carry out linearized 
fu simulations in which the gravitational potential, ^, is 
evaluated according to Eq. ([25|) . In the linearized treat- 
ment, the Compton wavelength is assumed to be fixed 
by the background field and thus chameleon effects are 
not present. Therefore, the difference between the full fu 
simulations and the linearized /r simulations are wholly 
due to the chameleon effects. To avoid confusion with 
linearization of the density field, we will call these runs 
the "no-chameleon" simulations. 

For each simulation box configuration, we run multiple 
simulations with different realizations of the initial power 
spectrum in order to reduce finite sample variance. The 
actual number of runs for each configurations are pri- 
marily constrained by computational resources and are 
summarized in Table HI To further reduce the sample 
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FIG. 4: Evolution of power spectrum deviation for f^Q = |10~^| 
for the full fji simulation and the no-chameleon simulation. The 
appearance of a chameleon at a < 0.5 causes a large fractional 
change in the deviations at earlier epochs. As the deviations grow, 
this offset remains as a smaller fraction of the total. 



TABLE L Simulation type and number 



i^box [h-^ Mpc) 





\fm\ 


256 


128 


64 


#of 


10-^ 


5 


5 


5 


boxes 


10-^ 


5 


5 


5 




10-^ 


5 


5 


5 




(GR) 


5 


5 


5 



Spatial Resolution {h~^ Mpc) 0.5 0.25 0.125 

kN/2 (h Mpc"^) 0.79 1.57 3.14 

Mass Resolution (lO^^/i"^ Mo) 27.6 3.45 0.431 



variance, we average the difference of the statistics per 
simulation from the ACDM run using the same realiza- 
tions of the initial conditions. 



B. PoAver spectrum results 

We start with the pure ACDM fno = simulations 
which serve as the baseline reference for comparison with 
the other cases. In Fig. [H we show the average power 
spectrum in the simulations of the three box sizes. The 
results from the various box sizes converge below about 
half the Nyquist wavenumber of the individual boxes. 
For comparison we also show the fit of Smith et al. [3y] , 
which scales the linear theory predictions (also shown) 
into the non-linear regime. The simulations converge to 



the accuracy of the fit again at about half the Nyquist 
wavenumber. 

All power spectra shown in this work, unless otherwise 
specified, are volume-weighted bootstrap-averaged from 
the multiple simulations of the same cosmological and 
f{R) parameters. In the averaging, the power spectra 
above the respective half-Nyquist modes are disregarded. 
Similarly, the error bars represent the volume weighted 
bootstrap error, in which the individual power spectrum 
data points are assumed to be uncorrelated. We use 2000 
bootstrap samples to compute the mean and the error. 
When we show differences of power spectra, the differ- 
encing is performed before averaging. 

In the bottom panel of Fig. [TJ we show the rela- 
tive difference between our simulations and the Smith 
et al. power spectra. As expected, the mean simula- 
tion power spectrum matches the Smith et al. results to 
~ 10% for all scales of interest. The error bars at the low 
k end of the simulation spectrum are large due to the 
small number of largest scale modes in Lbox = 256/i~^ 
Mpc box. In the opposite end, the errors are larger due 
to large variation between the different realizations of 
^box = 64/i~^ Mpc boxes. Note that bootstrap errors 
are only a rough estimate of the true uncertainties in 
these regions where the sample size is small. 

Fig. [2] shows the power spectrum enhancement of the 
|/i?o| = 10~^ and |/i?o| = 10~^ runs relative to the 
ACDM runs. Note that the enhancement of the power 
spectrum of the lensing potential <l> — ^ is identical by 
virtue of Eq. ([18]) . For comparison we also plot the linear 
theory predictions on the relative enhancement [1^, [3^ 
and the no-chameleon Jr simulation. In both the linear 
theory and the no-chameleon results, the force modifica- 
tion is completely determined by the background field. 
Thus on all scales smaller than the Compton wavelength 
in the background there is a sharp enhancement of power. 
Linear theory tends to overestimate the enhancement due 
to its neglect of mode coupling which makes the power 
at high k dependent on scales out to the non-linear scale 
where the effects are weaker. 

Since non-linear structures can only make the Comp- 
ton wavelength smaller, the no-chameleon and linear pre- 
dictions show the maximum scale out to which there are 
deviations from ACDM. However in the full Jr simula- 
tion the chameleon can dramatically change the results at 
high k if the background field amplitude is small enough 
to be overcome by the gravitational potentials of col- 
lapsed objects as discussed in ^ 

For the |/i?o| = 10~^ case the enhancement is reduced 
by a factor of 4 over the no-chameleon simulations at 
k ^ Ih Mpc~^. Previous simulations of modified force 
laws have all been in models where there is no chameleon 
and no way to hide deviations from ordinary gravity from 
local measurements [20|, [III, [H, [H, [^ • Note that even 
these reduced ~ 1 — 10% enhancements of power are po- 
tentially observable in next generation weak lensing sur- 
veys {e.g. [3^, [40]). As expected from the discussion in 
^n\ for the |/i^o| = 10~^ model typical gravitational po- 
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FIG. 5: Same as in Fig. [3] but at a = 0.3 for |/i?o| = 10~^ with the field values shown as a fraction of the background field at that 
time fpi{a = 0.3). The field amplitude is suppressed in deep potential wells exhibiting a chameleon that suppresses the early growth of 
structure. The shading scale range for the various fields in this case is also identical to that of Fig. [3] 
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FIG. 6: Power spectrum deviations from ACDM of the full /r 
simulations vs. the scaling predictions of Smith et al. employing 
linear fji calculations. Note that in all cases the scaling predictions 
fail to capture the deviations at high k due both to the change in the 
abundance and profiles of collapsed objects in the fji simulations 
and the chameleon mechanism. 



tentials of order 10~^ — 10~^ cannot overcome the back- 
ground field and the chameleon impact is greatly reduced. 
By taking slices through the simulations we can see the 
effects of the environment dependence of the chameleon 
(see Fig. [3]). In the |/i^o| = 10"^ run, the fn field is 
suppressed in the high density regions which correspond 
to deep potential wells. Notice that in low density re- 
gions the force law is still modified and this accounts for 
the small enhancement of power over the ACDM case 
that persists to high k in Fig. [2l Thus two identical 



sets of objects separated by the same distance will feel 
different forces depending on whether they are located 
in an overdense or underdense region. Generic tests of 
gravity such as the comparison between dynamical and 
lensing mass are predicted to produce null results in suffi- 
ciently overdense regions despite the ~ 1 — 10% enhance- 
ment of power shown in Fig. [2l Conversely, even with a 
chameleon, substantial modifications to the gravitational 
force law can appear in voids. 

In the |/i^o| = 10~^ run, the field remains stiff and 
at its background value across the whole volume leading 
to changes in the force law everywhere today. Even in 
this run, the power spectrum is still suppressed compared 
with the no-chameleon case (see Fig. [2]). This suppression 
comes about because the background field value was sub- 
stantially smaller at high redshift. Structures that form 
during these epochs are again affected by the chameleon 
and leave an impact at 2: = as hierarchical structure 
formation progresses. This can be seen in the evolu- 
tion of the power spectrum deviations in Fig. [H The 
impact of the chameleon at high z is fractionally large 
where the overall deviation is small. The offset remains 
roughly constant at more recent epochs as the overall 
deviation increases. In Fig. [5] we show slices through 
the simulation at a = 0.3 that reveal the presence of 
the chameleon in deep potential wells at that time (c/. 
Fig. [3]). The suppression of power spectrum enhance- 
ment at k > 0.7/1 Mpc~^ results from the shift of 1-halo 
contribution to larger scales and the relative flattening of 
A^(/c), as shown in the halo model inspired treatment in 
[9]. 

Finally, we can assess how well the Smith et al [S^ 
scaling works in the full fji simulations. In Fig. [6] we 
show that this prescription fails to capture the deviation 
from ACDM at high k. This disagreement appears for 
all of the fno values and arises both from changes in the 
contribution of collapsed objects to the power spectrum 
and the presence of the chameleon effect. 



In fact, the whole concept of the hnear power spectrum 
determining the non-hnear power spectrum at the same 
epoch that is shared by Smith et a/., the halo model, 
and other linear to non-linear scaling relations, is flawed 
in the context of a modification to gravity that evolves 
with redshift. We have seen that the precise form of the 
z = power spectrum in the |/i^o| = 10~^ runs depends 
on the presence or absence of a chameleon at a higher 
redshift. This information is not directly encoded in the 
linear power spectrum. 



IV. DISCUSSION 

We have carried out the first cosmological simulations 
of f{R) models for cosmic acceleration that exhibit the 
chameleon mechanism. The chameleon mechanism in- 
volves a non-linear field equation for the scalar degree 
of freedom that suppresses the range of the gravitational 
force modification or Compton scale in the deep grav- 
itational potential wells of cosmological and astrophys- 
ical structure. We have here focused on its impact on 
the matter power spectrum or equivalently the potential 
power spectrum relevant for weak lensing surveys. While 
we have simulated only one particular functional form of 
/(i?), we expect that the qualitative behavior of other 
models that exhibit a chameleon behavior to yield simi- 
lar results once scaled to the appropriate Compton scales 
and field amplitudes. 

In the absence of the chameleon mechanism, gravita- 
tional interactions would have an enhancement of a fac- 
tor of 4/3 on all scales smaller than the Compton scale in 
the cosmological background eventually leading to order 
unity enhancements in the power at high wavenumber. 
The chameleon mechanism turns on when the depth of 
the local gravitational potential becomes comparable to 
the field amplitude in the background. We have shown 
through otherwise identical simulations of structure with 
the chameleon mechanism artificially turned off that once 
the chameleon appears, it causes a substantial reduction 
of the enhanced power. For example, for a field ampli- 
tude of |/i^o| = 10~^ the change in the enhancement of 
power at /c ~ lh~^Mpc is a factor of four. 



Even in cases where current cosmological structures 
do not possess a chameleon {\fRo\ ^ 10~^), there still is 
an impact on the power spectrum due to evolutionary ef- 
fects. In f{R) models where the field amplitude decreases 
with curvature, the chameleon can appear again at high 
redshift when the building blocks of current structure 
were assembled. 

Scaling relations which take the linear power spectrum 
and map it into the non-linear regime qualitatively mis- 
estimate the non-linear power spectrum in several ways. 
For example, the Smith et al. |36] prescription assumes 
that the non-linear power spectrum depends only on 
the shape of the linear power spectrum near the non- 
linear scale. This prescription fails to describe both the 
chameleon mechanism and the change in the structure 
and abundance of collapsed objects leading to a severe 
misestimate at high k. 

A halo model can in principle do better to model these 
effects but simple prescriptions that scale the mass func- 
tion to the linear variance and leave halo profiles un- 
changed also do not describe the non-linear effects to 
sufficient accuracy ( cf. [9|] ) . In the next paper of this se- 
ries, we intend to study the impact of f{R) modifications 
on halo properties. 
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